Mixed Models II

Eva Freyhult

NBIS, SciLifeLab

June 11, 2025

Nested and crossed designs

Grouping structures can be complex and consist of more than one level of grouping.

Nested: lower level units are only associated with one higher level unit

Crossed: lower level units are associated with multiple higher level units

Nested designs

In many cases, we have a hierarchical structure in our data, where observations are nested within groups. For example, patients within clinics and clinics within regions.

graph TD
    A[Region 1] 
    A --> B[Clinic 1]
    A --> C[Clinic 2]
    A --> D[Clinic 3]
    
    B -.-> M1[Patient 1]
    B -.-> N2[Patient 2]
    B -.-> O3[Patient 3]

    C -.-> P4[Patient 4]
    C -.-> Q5[Patient 5]
    C -.-> R6[Patient 6]
    
    D -.-> S7[Patient 7]
    D -.-> T8[Patient 8]
    D -.-> U9[Patient 9]
    
    classDef region fill:#f9f,stroke:#333,stroke-width:2px;
    classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
    classDef patient fill:#cff,stroke:#333,stroke-width:1px;
    
    class A,E,I region;
    class B,C,D,F,G,H,J,K,L clinic;
    class M1,N2,O3,P4,Q5,R6,S7,T8,U9,V10,W11,X12,Y13,Z14,A15,B16,C17,D18,E19,F20,G21,H22,I23,J24,K25,L26,M27 patient;
    
    linkStyle default stroke-width:4px;

graph TD
    E[Region 2]
    E --> F[Clinic 4]
    E --> G[Clinic 5]
    E --> H[Clinic 6]
    
    F -.-> V10[Patient 10]
    F -.-> W11[Patient 11]
    F -.-> X12[Patient 12]

    G -.-> Y13[Patient 13]
    G -.-> Z14[Patient 14]
    G -.-> A15[Patient 15]
    
    H -.-> B16[Patient 16]
    H -.-> C17[Patient 17]
    H -.-> D18[Patient 18]
    
    I[Region 3]
    I --> J[Clinic 7]
    I --> K[Clinic 8]
    I --> L[Clinic 9]
    
    J -.-> E19[Patient 19]
    J -.-> F20[Patient 20]
    J -.-> G21[Patient 21]

    K -.-> H22[Patient 22]
    K -.-> I23[Patient 23]
    K -.-> J24[Patient 24]
    
    L -.-> K25[Patient 25]
    L -.-> L26[Patient 26]
    L -.-> M27[Patient 27]

    classDef region fill:#f9f,stroke:#333,stroke-width:2px;
    classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
    classDef patient fill:#cff,stroke:#333,stroke-width:1px;
    
    class A,E,I region;
    class B,C,D,F,G,H,J,K,L clinic;
    class M1,N2,O3,P4,Q5,R6,S7,T8,U9,V10,W11,X12,Y13,Z14,A15,B16,C17,D18,E19,F20,G21,H22,I23,J24,K25,L26,M27 patient;
    
    linkStyle default stroke-width:4px;

Nested model

By using a nested model, we can account for the hierarchical structure of the data and identify the sources of variability in our data.

Nested models in R

In R and lme4, random effects for nested groupings can be included. Random intercepts for both regions and clinics nested within regions can be included as follows;

nested_model <- lmer(y ~ x + (1 | region/clinic), 
                     data = data_int_nested)

(1 | region/clinic) specifies that we have random intercepts for clinics nested within regions

Note, (1|region/clinic) is identical to writing (1|region) + (1|region:clinic).

A third way to implement the same model is to make the nesting implicit by coding the lower levels (here clinics) uniquely, so that a clinic name only exist in one region. With such a naming scheme, we can use the notation (1|region) + (1|clinic).

Nested models in R

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | region/clinic)
   Data: data_int_nested

REML criterion at convergence: 1689.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.57686 -0.54281 -0.05204  0.44631  3.07737 

Random effects:
 Groups        Name        Variance Std.Dev.
 clinic:region (Intercept)  473.98  21.771  
 region        (Intercept) 2765.06  52.584  
 Residual                    50.84   7.131  
Number of obs: 240, groups:  clinic:region, 12; region, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept) 105.7145    27.2460   3.880
x             0.3223     0.1617   1.993

Correlation of Fixed Effects:
  (Intr)
x -0.124

Crossed designs

In a crossed design, an observation is associated with multiple groups. For example, patients may be treated by different doctors at several clinics.

graph TD

    A[Clinic 1]
    E[Clinic 2]
    I[Clinic 3]

    B[Doctor 1]
    C[Doctor 2]
    D[Doctor 3]

    M1[Patient 1]
    M2[Patient 2]
    M3[Patient 3]
    M4[Patient 4]
    M5[Patient 5]
 

    A -.-> B
    A -.-> C
    A -.-> D
    E -.-> B
    E -.-> C
    E -.-> D
    I -.-> B
    I -.-> C
    I -.-> D
   
    %% Patients connected to clinics to show they can come from any clinic, implying crossed with regions
    B --> M1
    B --> M2
    B --> M3
    B --> M4
    B --> M5
    
    C --> M1
    C --> M2
    C --> M3
    C --> M4
    C --> M5
    
    D --> M1
    D --> M2
    D --> M3
    D --> M4
    D --> M5
    
    %% Additional patient connections omitted for brevity

    classDef region fill:#f9f,stroke:#333,stroke-width:2px;
    classDef clinic fill:#ccf,stroke:#333,stroke-width:2px;
    classDef patient fill:#cff,stroke:#333,stroke-width:1px;
    
    class A,E,I region;
    class B,C,D,F,G,H,J,K,L clinic;
    class M1,M2,M3,M4,M5 patient;
    
    linkStyle default stroke-width:4px;

Crossed models

A crossed design is coded like this in R:

y ~ x + (x | clinic) + (x | doctor)

Model selection in mixed models

  • Mixed models can become complex and model selection can be used to identify the best model.

  • Likelihood ratio tests (LRT) are often used to compare nested models, where one model is a simpler version of another. LRT can be used also for mixed models.

  • Mixed model have both random and fixed effects variables. We can use LRT to compare models with different fixed effects structures, but also different random effects structures.

What is LRT?

The likelihood is a measure of how well a model explains the data.

The likelihood ratio test compares the likelihoods of two models: the null model, which is a simpler model, and the alternative, more complex, model.

\[\lambda = - 2\ln{ \frac{L(restricted\,model)}{L(full\, model)}} \] The test statistic follows a \(\chi^2\)-distribution with degrees of freedom equal to the difference in the number of parameters between the two models.

A significant difference in likelihoods suggests that the more complex model explains the data better than the simpler model.

For mixed models, we can use the anova() function to perform LRT. We will start by testing the random effects structure, without change of fixed effects. Then we will test the fixed effects structure, keeping the random effects fixed.

ML and REML

In mixed models, we can fit mixed models using two different methods, Maximum Likelihood (ML) or Restricted Maximum Likelihood (REML).

ML has the advantage that it can be used to compare models with different fixed effects structures, but it gives biased estimates of the variance components. More about this in Payam’s maths section.

In summary, use

  • ML when comparing models with different fixed effects structures (set refit = TRUE in anova())
  • REML when comparing diferent random effects in models models with the same fixed effects structure (set refit = FALSE in anova()). Also, always use REML for the final model.

Random Effects Comparison

Sleep study data, random intercept.

sleepall <- sleepstudy |> filter(Days>=2) |> mutate(Days=Days-2)
m1 <- lmer(Reaction ~ Days + (1 | Subject), data = sleepall)
summary(m1)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
   Data: sleepall

REML criterion at convergence: 1430

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6261 -0.4450  0.0474  0.5199  4.1378 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1746.9   41.80   
 Residual              913.1   30.22   
Number of obs: 144, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  267.967     10.871   24.65
Days          11.435      1.099   10.40

Correlation of Fixed Effects:
     (Intr)
Days -0.354

Random Effects Comparison

Sleep study data, random intercept and random slope for Days.

m2 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepall)
summary(m2)
Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepall

REML criterion at convergence: 1404.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0157 -0.3541  0.0069  0.4681  5.0732 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 958.35   30.957       
          Days         45.78    6.766   0.18
 Residual             651.60   25.526       
Number of obs: 144, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  267.967      8.266  32.418
Days          11.435      1.845   6.197

Correlation of Fixed Effects:
     (Intr)
Days -0.062

Random Effects Comparison

Compare the two models using a likelihood ratio test.

anova(m1, m2, refit = FALSE)
Data: sleepall
Models:
m1: Reaction ~ Days + (1 | Subject)
m2: Reaction ~ Days + (1 + Days | Subject)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
m1    4 1438.0 1449.9 -715.01   1430.0                         
m2    6 1416.1 1433.9 -702.05   1404.1 25.926  2  2.346e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fixed Effects Comparison

We can also compare models with different fixed effects structures. For example, we can compare a model with a fixed effect of Days to a model with a fixed effect of Days and Sex.

m2 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepall)
m3 <- lmer(Reaction ~ Days + sex + (1 + Days | Subject), data = sleepall)
anova(m2, m3, refit=TRUE)
Data: sleepall
Models:
m2: Reaction ~ Days + (1 + Days | Subject)
m3: Reaction ~ Days + sex + (1 + Days | Subject)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
m2    6 1425.2 1443.0 -706.58   1413.2                       
m3    7 1422.3 1443.1 -704.14   1408.3 4.8826  1    0.02713 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

P-value calculations

As you have noticed lme4 does not report p-values for the fixed effects.

In linear regression, p-values are calculated based on the t-distribution of the estimated coefficients. The calculation of denomicator degrees of freedom is straightforward since all observations are independent.

In mixed models, the situation is more complex due to the grouping structure. It is not straightforward to calculate the degrees of freedom for the fixed effects, since the observations are not independent.

There are however some approximation methods available to calculate the degrees of freedom and hence p-values for the fixed effects in mixed models.

We can use the lmerTest package, which extends lme4 to include p-values for fixed effects, based on e.g. Satterthwaite’s method.

lmerTest

library(lmerTest)
mm <- lmer(conc ~ age + (1 | group), data = df)
summary(mm)
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: conc ~ age + (1 | group)
   Data: df

REML criterion at convergence: 304.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.75481 -0.89181  0.09125  0.70848  1.85339 

Random effects:
 Groups   Name        Variance Std.Dev.
 group    (Intercept) 2224.44  47.164  
 Residual               85.47   9.245  
Number of obs: 40, groups:  group, 4

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 241.3334    26.0615   4.2275   9.260  0.00058 ***
age          -1.8293     0.3014  36.1653  -6.068 5.54e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
    (Intr)
age -0.422
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy

REML criterion at convergence: 1743.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4634  0.0231  0.4634  5.1793 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.10   24.741       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  251.405      6.825  17.000  36.838  < 2e-16 ***
Days          10.467      1.546  17.000   6.771 3.26e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
Days -0.138

Predicting random effects